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Chapter 1 
Introduction 



From a mathematical point of view, quantum mechanics can be regarded in many respects 
as an eigenvalue problem. Typically, one has to to calculate eigenvalue and eigenvectors 
of Hamiltonians. Very often, especially in low temperature problems, the knowledge of the 
ground state and the first few excited states yields considerable insight into the physics of a 
given system. Symmetry properties and the known quantum numbers of the desired state can 
be used to reduce the Hilbert space. However, even a large reduction factor will eventually be 
overcome by the exponential growth. This means that requirements for memory and running 
time of relatively small systems can be prohibitive, even for the best computers, when a full 
diagonalization of a finite cluster is attempted. Typically, many degrees of freedom have 
to be integrated out of the original problem to make it accessible to present day computer 
capacities. 

This notes consist of an introductory course to the Lanczos Method and Density Matrix 
Renormalization Group Algorithms(DMRG), two among the leading numerical techniques 
applied in studies of low-dimcnsional quantum models. Another important numerical tech- 
nique, the Quantum Monte Carlo, is discussed by R. R. dos Santos in his notes for this 
School. The numerical approach allows for a direct and unbiased calculation of physical 
properties for finite clusters, from which a phase diagram can be constructed. 

For those outside the field of one- or two-dimensional quantum systems it may seem 
that working on a such low- dimensional is completely irrelevant to three-dimensional reality. 
However, there are in fact may reasons for working in lower-dimensional physics, in both 
statistical mechanics as well as other fields. 

First, restricting study to one dimension continues to demonstrate itself as an efficient 
and effective laboratory for the development and consideration of new theoretical ideas, 
many of which arc intended for application to real higher-dimensional. Both thermal and 
quantum fluctuation effects are larger in lower dimensions. 

Secondly, there are in fact real systems which demonstrate a high degree of low- dimensionality, 
at least in condensed matter physics. 

Thirdly, one-dimensional systems share the possibility with higher-dimensional systems 
of undergoing quantum phase transitions as some parameter is varied while holding the 
temperature fixed at T = 0. In experimental situations, this parameter could be pressure or 
doping, for example, such as in the case of high temperature superconductors where, upon 
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the variation of doping, the ground state ranges from antiferromagnetic to superconducting. 
This manuscript is organized in the following way. In Section 2 the Hamiltonians for the 

Hubbard, t — J and Heisenberg models are defined. This models and their generalizations 
are aimed to describe the microscopic features of matter. Besides that, they present rich 
matematical aspects. In Section 3, the idea of studying the models on clusters of a finite size 
in order to extract their physical properties is discussed. The important role played by the 
model symmetries is also examined. In Section 4 the Lanczos method is presented and some 
of its most popular generalizations discussed. In Section 5 we start by examining the basic 
ingredients of the DMRG method and presenting the two kinds of DMRG algorithms. Then, 
the way measurements of observables are performed in DMRG is discussed and, finally, the 
most important additions to the original method are briefly mentioned. Section 6 contains 
some final comments and the acknowledgments. 
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Chapter 2 
Models 



In this Section we present some representative quantum models. These models and their 
generalizations have wide applicability in statistical mechanics as well as in condensed matter 
physics. 



A starting point for the microscopic approach to electron motion in crystals may be obtained 
by analyzing the energy levels of the atoms involved. Lattice models are used, because the 
atomic structure in the physical systems determines the possible places at which the electron 
can be found. If two neighboring atoms have overlapping orbitals with very similar energies, 
the orbitals can hybridize and allow electrons to travel from one atom to the other. On the 
other hand, the repulsion between electrons is very strong due to their charge. The simplest 
approximation for the interaction between the electrons is to restrict it to the case when both 
electrons are on the same site (or atom). On-site Coulomb repulsion and nearest- neighbour 
(n.n.) hopping are already the terms of the Hubbard model[4, 5], which can be characterized 
by the Hamiltonian 



where (cjo-) is a creation (annihilation) operator for an electron of spin o" =1, | in a 
Wannier orbital at lattice site i and (ij) denotes n.n. pairs. In the sum, bonds (ij) are 
summed over only once each. Here t is the matrix element for tunneling from one lattice site 
to a neighbouring one, i.e. the overlap of n.n. electron wavefuctions. The letters h.c. denote 
the Hermitian conjugate of the immediately preceding term. The fermionic operators obey 
the anticommutation relations 



2.1 The Hubbard Model 




(2.1) 




(2.2) 



and 



{Cia, Cjfi] — 0. 



(2.3) 
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Also, riifj — cj^Cifj is the number operator for electrons of spin a at site i, and U the Coulomb 
repulsion. The charge or number operator at site i is 

"-i = = + nil- (2.4) 

The states of the model are given by specifying the four possible configurations of each site 
(its level can either be empty, contain one electron with either of two spins, or two electrons 
of opposite spins) on a lattice made of L sites. 

However, the knowledge of the parameters t and IJ is not enough to characterize the 
system. One also needs to know: 

• The dimensionality of the system, i.e., whether it is a one- dimensional (ID) chain, a 
two-dimensional (2D) plane or a full three-dimensional (3D) system. 

• The geometry of the lattice. The lattice structure can introduce frustation, and its 
symmetry properties greatly influence the behaviour of electrons traveling them. 

• The boundary conditions (BC), i.e., whether the lattice is opened, closed or has some 
special condition at its surface. 

• The filling or density, 

n=\Y.M. (2.5) 

where (...) denotes a ground state expectation value. 

• The temperature T which is the third energy scale along with t and IJ . 

The first term in Hamiltonian Eq. (2.1) is called hopping term and the second one 
Coulomb term. The hopping term alone can be shown to lead to a conventional band 
spectrum and one-electron Bloch levels in which each electron is distributed throughout the 
entire crystal (a metal). The Coulomb term alone would favor local magnetic moments, 
since it suppress the possibility of a second electron (with oppositely directed spin) at singly 
occupied sites (an insulator). When both terms are present the competition among them 
brings about a transition between the metallic phase and the Mott insulating phase [6]. 

2.2 The t - J Model 

If the Hubbard model is considered in the limit where U/t\s large, the strong coupling limit, 
the number of doubly occupied sites is small. This leads to the derivation of effective models, 
most prominently the t — J model[7]. In the t — J model, the Hubbarb model complexity is 
reduced by projecting out the states with double occupancy. Using this procedure one gets 
the t — J Hamiltonian 

HtJ = Z] (claCja + ^-C.) + JYli^^i' 4 ~ ' (2-6) 

(y>o- {ij) 
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where 

Si^J2(^ial^apCip (2.7) 

a/3 

is the electron spin operator at site i, with 



a = [a ,CT ,CT 




(2.8) 



being the Pauh matrices. The tilde on the c-operators in the hopping term refers to the fact 
that a creation operator cannot introduce an electron on a site where another electron, even 
of the opposite spin, is already located. Formally, one would express this as 

4a = 4a(l - '^i-a) and, thus, Cia- = (1 - ni^-a)cia- (2.9) 

In the derivation of the t — J model from the Hubbard model, one takes in account 
intermediate double occupancies between states which are suppressed by their high energy 
but help to lower the kinetic energy by making the electrons more mobile. In such processes, 
an electron hops on an already occupied site, then cither electron can hop back to the original 
site. Double occupancies are forbidden in the t — J model, but the physics can be included by 
adding an additional term. Since the intermediate processes in the Hubbard model are only 
possible if the electrons on the adjacent sites have opposite spin, this term can be described 
by an interaction that favors the singlet state compared to the triplet state. This is precisely 
what the second term in Hamiltonian Eq. (2.6) docs. To show this, we note that Si ■ Sj has 
two eigenvalues: if the two electrons are in a singlet state, the eigenvalue is -3/4, and if the 
electrons form a triplet, the energy is 1/4. If the effect of the |njnj-term is included, one 
gets the following energies from the J-term: 

one or both sites unoccupied 
spins on both sites forming a triplet state 
spins on both sites forming a singlet state —J. 

Therefore, the term effectively lowers the energy for states in which two electrons with 
opposite spin are on adjacent sites. Prom second order perturbation theory of the Hubbard 
model we obtain J = 4:t'^/U[7]. 



2.3 The Heisenberg Model 

In the case of half-filling (n = 1), when Ui = 1 for all sites, hopping becomes impossible in 
the t — J model. In addition, the |njnj-term is reduced to a constant. Thus, both terms 
can be neglected and the t — J model at half fiUing is a spin 1/2 Heisenberg model. 

It is important to remember that this is an effective model, i.e., although the spins in the 

— * — * 

model are localized they are meant to describe a system of mobile electrons. The Si ■ Sj 
interaction is called a spin exchange interaction. 
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All these models were originally designed to describe physical systems. The Heisenberg 
model was introduced before the Hubbard model and is not restricted to just spin 1/2 
systems. The restriction in the above case resulted from the fact that electrons are spin 1/2 
fermions. 

A possible generalization for these models is to extend the n.n. summation to more 
distant neighbours. 
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Chapter 3 

The Study of Finite Size Systems 



Having set up Hamiltonians which are beheved to contain the physics we are interested 
on, we are left with the formidable task of calculating measurable quantities to gain some 
understanding of the models. 

Here we choose to work with two unbiased numerical techniques namely, Lanczos and 
DMRG. These methods are unbiased {i.e., in contrast to mean-field based approaches) in the 
sense that they do not make an initial assumption on the nature of the ground state of the 
system. However, as we will see later, each of these numerical methods has its limitations. 

In studying finite size systems the general idea is to construct a matrix representation 
of the Hamiltonian for a given number of sites (system size). The Hamiltonian is then 
diagonalized in order to obtain the spectra and calculate measurable quantities, e.g., spin and 
charge correlation functions. We repeat that for systems of different sizes and extrapolate the 
results toward an infinite size system {i.e., towards the thermodynamic limit). This agenda 
is usually well succeeded provided that we have enough data for a reliable extrapolation. The 
number of lattices with different sizes needed for the extrapolation procedure to converge 
depends heavily on the model being studied and even on the set of parameters being used for 
a given model. The Finite Size Scaling theory[8] is specially useful when critical behaviour 
is present since, in this case, a strong dependence of the physical quantities on the system 
size is expected[9]. However, in many cases a fairly good ideia of the properties in the 
thermodynamic limit can be achieved by examining just a few system sizes. The models we 
are dealing with involve many parameters {e.g., n, t, U, and J in the models of Section 2). 
In order to build up a phase diagram we have to systematically cover the model parameter 
space to see how the physical quantities depend on each of these parameters. 

We can roughly divide the above task in two steps: 

• Build up a representation for the Hamiltonian, diagonalize it, and calculate the mea- 
surable relevant quantities. 

• Interpret the results and construct the phase diagram. 

The second step strongly depends on the model being studied and on what kind of physics 
we are looking at. If the subject being studied is an advanced topic in physics (as it usually 
is), then a considerable experience in that field of research might be necessary to carry that 
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step out. On the other hand, the first one is, in principle, much more simple and involves the 
knowledge of basic concepts in quantum mechanics and numerical analysis; topics covered 
in most undergraduation courses in physics. The Lanczos and DMRG are tools for carrying 
out the first step. To understand why such special techniques are needed let us suppose 
the Heisenberg model is to be studied on a lattice of L sites. Each site has two possible 
states: spin up and down. A lattice with L sites has 2^ states and this is the dimension of 
the Hamiltonian matrix. Similarly, for the t — J and Hubbard models we have 3^ and 4^, 
respectively. Due to this exponential growth with L even small lattices (typically 10 sites or 
so) generate Hamiltonians too big to be handled by present-day computers using standard 
diagonalization algorithms. 



3.1 Symmetries 



In the process of constructing a representation for the Hamiltonian it is very useful to take 
advantage of the model symmetries. Many models, including those presented in Section 2, 
exhibit conservation of total spin number, total spin in the 2;-direction and total charge, i.e., 



H,{SY =[H,S^] = [H,N]^0, 



where H is the model Hamiltonian and 



i 



(3.1) 

(3.2) 
(3.3) 



In addition, these operators also commute with one another, i.e., 



= 0, 



(3.4) 



so that the eigenvalues of H, S, S^, and are simultaneous good quantum numbers, which 
we will denote simply by E, S{S + 1), S^, and N respectively.^ In the numerical treatment 
of a given model it is possible to consider eigenstates which simultaneously diagonalize 
H and all operators associated to its symmetries. We do this by choosing to work in a 
representation in which the symmetry operators are always diagonal, selecting a subspace or 
sector of Hilbert space with particular eigenvalues of those operators, and diagonalizing H 
in this particular sector. As exemplified below, total spin in the 2;-dircction and total charge 
are easily implemented. However, total spin quantum number (S)'^ is much more difficult to 
specify, and also quite difficult to measure, since in terms of the fundamental operators Qo- 
it is extremely non-local: 



(3.5) 



^It will be clear from the context whether by or N we mean the operator or its eigenvalue. 
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where i and j both run throughout the entire lattice, i.e., they do not just represent n.n. 
terms. However, in cases in which choosing the direction of quantization in the z-axis is 
arbitrary, i.e., in which there is full rotational (or SU{2)) symmetry, the eigenvalues of 
H, N, and S"^ will be independent of the eigenvalues of S^, which may assume any value 
ranging from —S to S. This can be quite useful in determining the total spin quantum 
number S of the ground state if it is only possible to specify the projection in the method 
of investigation. 

Let us denote by E{S^) the ground state in the sector specified by the spin projection 
S'. If E{S') < E{S' + 1) and E{S') = E{S' - 1) = ... = E{S^J, where 

_ / if -S^ integer , . 

^min I i jf half-odd-integer, ^ ^ 

then the absolute ground state, (i.e., that in the Hilbert space unrestricted by the specifica- 
tion of S^) contains a state of spin and no states of any higher spin. 

Strongly correlated electron models also often exhibit a particle-hole symmetry. This 
symmetry relates the creation of an electron to its destruction in the following way. Consider 
the transformation 

PH : ^ (-l)VL. (3.7) 
Under this transformation, the n.n. hopping terms transform according to 

PH : cj^Ci+i^a + h.c. -Ciacl^-^^^ + h.c. = cl^i^^Cia + h.c. = cl„Ci+i^a + h.c, (3.8) 

remaining unchanged under this transformation. However, the number operators transform 
according to 

PH . Tliu — C^j^^CiQ- > Ci(jC^^ — 1 C^^Cifj — 1 Tlicr (3-9) 

FR-.rii 2-ni (3.10) 
FR : N 2L-N, (3.11) 

and similarly the conduction electron spin operator transforms according to 

Fli:s, = {slslst) ^ (-sls!,-st) = m (3.12) 
PH : ^= (5^S^5") ^ {-S\Sy,-S') ^RS (3.13) 
PR : S\ (3.14) 

Here R represents a spin rotation by tt about the y-axis. If this particle-hole symmetry 
applies to the full Hamiltonian, i.e. not only to its non- interacting kinetic part as shown 
here but also its interaction terms, then there will be a one-to-one correspondence between 
eigenstates with quantum number {E, S, , N) and {E, S, —S^ ,2L — N). In particular, if 
we wish to determine the properties of a system away from half-filling, we need only do this 
below half-filling and the physics above half-filling will be identical. 

In order to reduce finite size effects it is common to work with periodic boundary con- 
ditions (PBC) to eliminate the lattice surface in given a direction. In this case, the system 
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gains translational invariance in that direction^. This symmetry express itself through the 
commutation of H with the translational operator T, which can be defined as 



T\did2...dL-idL) = \d2ds...dLdi), (3.15) 

where 

\did2...dL-idL) = \di) (8) \d2)^, OjciL-i) ^ W (3.16) 

and di is the state at site i {e.g., for the Hubbard model di can be either one electron with 
spin up or down, two electrons with opposite spins, or an empty site). In Eq. (3.15) we have 
assumed, for the sake of simplicity, a one dimensional lattice of size L. In higher dimensions 
translational invariance in each direction can be treated separately. The eigenvalues of T 
are e^^^, A; = 0, 1, 2, L — 1, which we will label by the number k. An eigenstate of T 
with eigenvalue k is given by 



C 



l + U + + u'^ + ... + u^-^ 



\did2...dL-idL), (3-17) 



where U — e'"r'=T and C a normahzation constant. If 0^ = then is not possible to construct 
a state with momentum k from |did2---c^L-iC^L)- Note that S, N, S^, and k are simultaneous 
good quantum numbers. 

Another symmetry shared by all the models introduced and Section 2 and many others 
is the charge-conjugation symmetry. This symmetry implies a one to one correspondence 
between eigenstates with quantum number {E, S, S^, N) and {E, S, —S^, N) and is present 
if the Hamiltonian commutes with the parity operator, whose effect is to flip the particle 
spin (|TiOOT)~^liTOO|)). In the presence of charge-conjugation symmetry we can 
choose to work only with 5*^ < or only with > 0. 

Let us consider the t — J model on a chain with four sites under PBC. Each site can 
be in one of the following states: one electron with spin up (| t)) or down (| J,)), or empty 
(|0)). The dimension of total Hilbert space is 3^ = 3^ = 81. We divide the Hilbert space 
in sectors labeled by the quantum numbers A^, S^, and k. Below, we denote each sector by 
[N, S^, k, (dimension of the sector)] and write down the states for some illustrative cases. 



[0,0,0, 1] 

[1,0.5,0,1] [1,0.5,1,1] [1,0.5,2,1] [1,0.5,3,1] 
[2,0,0,3] 



4ao,3] = ^(|0 0iT) + |0iT0) + |iT0 0) + |T0 0i)) 
4ao,3] = ^(|OTOi) + |toio) + |oiOT) + |iOTO)) 



In fact, a cyclic BC is enough to gain translational invariance. 
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[2,0,1,3] 



[2,0,2,3] 



[2,0,3,3] 



[2,1,0,2] 



[2,1,1,1] 



[2,1,2,2] 



[2,1,3, 1] 



4ai,3] = ^(|ooTi)+i|OTio)-|Tioo)-i|iooT)) 
4ai,3] = ^(|oon)+^|oiTO)-UTOO)-^|TOOi)) 
4ai,3] = ^(|OTOi) + i|TOiO)-|oiOT)-iUOTO)) 



4a2,3] = ^(|oon)-|OTio) + |Tioo)-|iooT)) 

(.(2) 1 



4a2,3] = ^(|otoi)-|TOio) + |oiot)-UOTO)) 



.(1) 1 



c{2) 

[2,0,3,3] 2 



i'U] = ^(|oioi) + uoio)) 



\|0 0iT)-^|0iT0)-UT0 0) + ^|i0 0T)) 



0[2,O,3,3] = 2^'°^°^^"''^°^°^"'°^°^^+''^°^°^^ 



[iW] = 2(|0 0^^) + |0^^0) + l^^0 0) + l^0 0^)) 



1 



4'm,i] = i i) + i|0 i i 0) - u i 0) - i )) 



2 



(.(1) 1 



4l2,2] = ^(lOiOi)-liOiO)) 



[iW] = odO i i) -^|0 i i 0) - li i 0) +zU i )) 



2 
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[3,0.5,0,3] [3,0.5,1,3] [3,0.5,2,3] [3,0.5,3,3] 



. [3,1.5,0,1] [3,1.5,1,1] [3,1.5,2,1] [3,1.5,3,1] 
. [4,0,0,2] [4,0,1,1] [4,0,2,2] [4,0,3,1] 
. [4,1,0,1] [4,1,1,1] [4,1,2,1] [4,1,3,1] 

• [4,2,0,1] 

The example above give an ideia of how helpful symmetry implementation can be. Indeed, 
a further reduction in the Hilbcrt space can be achieved in some cases by considering a lattice 
reflection symmetry, namely [W, Htj] — 0, where 

W\did2...dL-idL) = \dLdL-i...d2di). (3.18) 

The eigenvalues of W are: +1 and -1. For instance, the sector [2,0,0,3] can be broken in 
two subspaces, indexed by = — 1 ([2, 0, 0, —1, 1]) and +1 ([2, 0, 0, 1, 2]), and given by: 

• [2,0,0,-1,1] 

• ([2,0,0,1,2]) 

<^[2A0,l,2] - ^l</^[2,0,0,3] +'/^[2,0A3]J 

^(2) _ J,(3) 

V^[2,0,0,l,2] — V'[2,0,0,3]- 

In higher dimensions, besides reflection with respect to different axes, rotations by several 
angles arc also available. Note that these operations might or might not commute with the 
translational operator depending on the value of k. For instance, reflection symmetry can 
not be employed to break sector [2, 0, 1, 3] in sub-sectors. 

By implementing the model symmetries we can significantly reduce the computational 
effort required for diagonalizing the Hamiltonian. In some cases we can study lattices large 
enough to reveal bulk properties using standard routines to diagonalizc each sector of the 
Hamiltonian. For instance, a 4x4 Heisenberg system can be fully diagonalized if most sym- 
metries discussed above are implemented[10]. Since this approach yields the full spectrum we 
can construct the exact partition function for the system and, therefore, obtain the thermal 
behaviour exactly at arbitrary temperature.^ 



^As discussed below, this is a commodity we will not have when working with Lanczos or DMRG methods. 
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Chapter 4 
Lanczos Method 



In this Section, an algorithm is presented which allows us to determine numerically the 
ground state and some excited states for Hamiltonian operators on finite clusters. The basic 
idea of the Lanczos method[ll, 12] is that a special basis can be constructed where the 
Hamiltonian has a tridiagonal representation. Once in this form the ground state of the 
matrix can be found easily using standard library subroutines such as Numerical Recipes or 



The tridiagonal matrix is constructed iteratively. First, it is necessary to select an arbi- 
trary normalized vector I'^o) in the Hilbert space of the model being studied. The overlap 
between the actual ground state |^o)i ^^^d the initial state I'^o) should be nonzero. If no a 
priori information about the ground state is known, this requirement is usually easily satis- 
fied by selecting an initial state with randomly chosen coefficients in the working basis that is 
being used. If some other information of the ground state is known, like its total momentum 
and spin, then it is convenient to initiate the iterations with a state already belonging to 
the subspace having those quantum numbers (and still with random coefficients within this 
subspace). 

After \iJjq) is selected, we define a new vector by applying the Hamiltonian H to the initial 
state. Subtracting the projection over \ipo), we obtain 



that satisfies (V'olV'i) = 0. Now, we can construct a new state that is is orthogonal to the 
previous two as. 



It can be easily checked that (V'olV'2) = — 0. The procedure can be generalized by 

defining an orthogonal basis recursively as. 



IMSL. 




(4.1) 






(4.3) 



where i — 0,1, 2, and the coefficients are given by 




(4.4) 
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supplemented by 60 — 0, — 0. In this basis, it can be shown that the Hamiltonian 

matrix becomes, 



ao 


bi 








bi 


ai 


b2 








62 


a2 


bs 








bs 


as 



\ ; ; ; ; ; 

i.e., it is tridiagonal as expected. Once in this form the matrix can be diagonahzed using 
standard hbrary subroutines. However, note that to diagonahze completely the model being 
studied on a finite cluster a number of iterations equal to the size of the Hilbert space (or of 
the subspace under consideration) is needed. In practice this would demand a considerable 
amount of CPU time. However, one of the advantages of this technique is that accurate 
enough information about the ground state can be obtained after a small number of iterations 
(typically of the order of ~ 100 or less). The ideia behind Lanczos method is a systematic 
improvement of a given variational state that is used to represent the ground state [14] . The 
procedure just described assumes H being an hermitian matrix. If that is not the case then 
a generalized algorithm is nccdcd[15]. 

The eigenvalues of (4.5) steadily approach the lowest eigenvalues of H and its eigenstates 
are expanded in the Lanczos basis \ipi). Each state I'^j) is represented by a large set of 
coefficients, when it is itself expanded in the basis selected to carry out the problem {e.g, 
the basis used in example of Section 3.1). Thus, in practice, it is not convenient to store 
each one of the {ipi) vectors individually, since such a procedure would demand a memory 
requirement equal to the size of Hilbert space sector multiplied by the number of Lanczos 
steps. A simple solution to this problem consists of running the Lanczos subroutine twice. 
For instance, if j^^o) = J2ifi\'^i): then in the first run the coefficients fi are obtained, and 
in the second the vectors are systematically reconstructed one by one and used to build 
up |\E'o) in the original basis. 

A common difficulty with the Lanczos method is that finite precision arithmetic causes 
the to lose their orthogonality. A consequence of that is the appearance of spurious 
eigenvalues (and eigenstates) in the spectra. One fix is to repeatedly reorthogonalize, which is 
too costly since, in this case, all have to be stored. Another is to partially reorthogonalize, 
and a third option is to ignore the problem. We can choose the later if only the extremal 
(lowest or highest few) eigenvalues are needed, which is often the case. A final way of 
avoiding the problem is to stop the Lanczos whenever the problem starts to appear, calculate 
the ground state, and use it as an initial state \ipo) to restart the Lanczos procedure. This 
resets the basis but keeps the information from the previous Lanczos running. By pushing 
this ideia to its limit, we can perform always just two Lanczos steps {i.e., work just with \iIjq) 
and iV'i)), diagonahze a 2x2 matrix, and use its lowest eigenstate as a new \ipo). This is the 
so-called modified Lanczos[16, 17]. The modified Lanczos converge slowly than the original 
Lanczos but has the convenience of having the ground state always at hand. An even more 
pedestrian technique is the power method[19], which consists of successively applying the 
Hamiltonian to the initial state until all excited states are filtered out and only the ground 
state remains. This procedure is the slowest in speed of convergence, but in simple problems 
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is enough and easy to program. 

One of the greatest appeals of the Lanczos method is the possibihty of calculating dy- 
namical properties of a given Hamiltonian in finite clusters[18, 12]. The technique permits 
accurate calculations of energy and momentum dependent dynamic correlation functions 
which are observable in scattering experiments, such as Neutron Scattering (spin dynamics) 
and Photoemission Spectroscopy which measures the spectral function of the system [13]. 

Since Lanczos yields a number of excited states another interesting possibility is the 
calculation of finite-temperature quantities. Examples of successful attempts in this front 
are Refs. [20] and [21]. 

The main limitation of Lanczos technique is the size of the clusters that can be studied. 
Recently, attempts have been made to reach larger cluster. The basic ideia is the following. 
If is a complete basis we are working with, then the ground state can be formally 
represented as 

l*o) = E^^I'^^)- (4-6) 

i 

In general, all \(f)i) contribute significatively for the sum but, in some cases, it might happen 
that several states have very small weight gi. In fact, for most models studied a very small 
percentage of the states in the basis dominate the sum in Eq. (4.6). This suggests that 
useful results can still be obtained if a large part of the Hilbert space is simply discarded. 
The truncated Lanczos[22] implements this idea by systematically enlarging and reducing 
the working Hilbert space in a controlled way and by keeping a fixed number of states 
(typically a million or so) in the basis. The number of states kept represent a very small 
percentage of the total Hilbert space but the algorithm is designed to search for and keep 
the dominant states in sum Eq. (4.6). It is also important to choose a basis which uses 
as much information about the system as possible. For every physical systems we should use 
a different basis depending on physical insight we have about the system. For instance, if 
we know that there is a tendency to dimerization then it is convenient to construct a basis 
of dimcrs. A well chosen basis leads to a better representation of |^o) for ^ given number of 
states kept in the truncation procedure. 

The efficiency of the truncation technique depends heavily on the model being studied. 
In particular, it seems suitable for problems with gaps in the spectrum (like a spingap) [23] . 
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Chapter 5 



Density Matrix Renormalization 
Group Algorithms 

5.1 Overview 

The basic agenda to overcome the system size hmitations is to use a basis in which the 
ground state can be represented by only few basis states. In other words, a procedure must 
be found to identify or construct the important states and neglect or discard all others so 
that the piece of the Hilbert space one operates on remains small. The truncated Lanczos, 
briefly discussed in the previous Section, is a possible approach to this agenda. It has the 
advantage of working with a basis formed by states that have an intuitive meaning so that 
the results can be easily interpreted. In addition, dynamical information can be obtained 
without difficulty. The DMRG technique we are about to discuss is an alternative approach 
to that agenda. The innovation of the DMRG is that it does not hold on to a specific basis, 
but optimizes the basis it uses in the steps leading to the calculation of the ground state. 
A disadvantage of the DMRG method is that the basis states chosen by the algorithm are 
not intuitive, and the description of the state requires the measurement of observables. For 
the measurement process, one needs a representation of the operators in the current basis. 
Consequently each operator that needs to be measured must be stored, and every time the 
basis is changed all of them have to be transformed. This is expensive in time and memory. 
Another disadvantage is that dynamical information cannot be easily obtained. 

Historically DMRG has its roots in the renormalization group approach pioneered by 
Wilson[24]. Considering a real-space blocking version for lattice systems of Wilson's original 
approach, the basic idea is to start with a small system that can be handled exactly. The 
system size is then increased without increasing the size of the Hilbert space until the desired 
system size is reached. 

Increasing the system without increasing the Hilbert space is typically done in two steps: 

• The system size is increased, and as a consequence, the Hilbert space grows at the 
same time. 

• The Hilbert space is truncated to its original size keeping the system size constant. 
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To characterize such a renormalization procedure two basic questions have to be an- 
swered: 

• How is the enlargement done ? 

• Which criterion to apply in the second step to differentiate between the states we will 
keep from those we will discard ? 

In Wilson's like approach, we start with blocks of a certain (small) size. In the first step, 
two such blocks arc linked to form a block which is twice as large. The Hamiltonian of this 
larger block is then exactly diagonalized and its eigcnstates are used as basis states. The 
criterion for keeping states is their energy, and only eigenstates whose energy lies below a 
certain threshold are kept. The states which are kept characterize the new block that is 
again linked to an identical block, and the process is iterated. 

This approach proved to be very effective for the Kondo model[24]. However, for other 
strongly correlated systems like those in Section 2 it was not successful [25]. The main reason 
for this failure lies in choosing the block eigenstates as the states to be kept. Since the block 
was not previously connected to the rest of the system ( another identical block in the case 
above ) its eigenstates have inappropriate features at the block ends, making them a poor 
choice as a basis to represent the ground state of a larger system, formed by putting together 
two (or more) blocks. This problem was pinpointed by White and Noack in Ref. [1] and 
an attempt was made to fix it by combining eigenstates from several different blocks under 
various BC. Let us see how DMRG fix this problem. 

5.2 Enlargement and Reduction in the DMRG Proce- 
dure 

As mentioned in Section 5.1, the two most important characteristics of a renormalization 

procedure are the way the system grows and how the decision is made on which states are 
kept in the Hilbert space by the truncation step. In this Section, these elements of the 
DMRG procedure will be discussed. How these elements are used in the global DMRG 
algorithm will be the subject of the subsequent Sections. 

Figure 5.1 shows the most important structures used in the DMRG algorithm. The 
elementary unit is a site, and is described by the states di {i = 1,...,D), in which the 
site can be found^. A block B{l^m) consists of a number of sites / and its Hamiltonian Hb 
contains only terms involving the sites inside the block. To represent B{l,m) and Hb we 
associate to them a m-dimensional basis where m is in general smaller than the full Hilbert 
space of the block. The states in the basis are grouped in symmetry sectors labeled by a 
set of quantum numbers {e.g., and N), which makes Hb a block-diagonal matrix. We 
also store the matrix elements of Hb between these states. The block is grown by adding 
a site to it, and together they form the enlarged block B"^. If . . . \bm) and \di) . . . \dr)) 

^Here, the index i is not labeling a site in a lattice but the states accessible to a given site. For instance, 
D = 4 and 2 for the Hubbard and Heisenberg model, respectively. 
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site # 



block B(l,m) 




enlarged 
block ^ 



Figure 5.1: Representation of the basic elements of the DMRG algorithm. The thick line 
connecting the block and the site represents all the interaction terms among them present 
in the Hamiltonian. 



represent, respectively, the basis of block and a site then the basis of the enlarged block can 
be constructed from the direct product 

\bl)^\b,)®\dj). (5.1) 

The dimension of the Hilbert space for B'^ is the product of the dimensions of the Hilbert 
space of B{1, m) and a site, i.e, m x D. A possible mapping of i and j-values onto /c-values 
in Eq. (5.1) is k = {i - 1)D + j. 









B(l,m) 




B'(l',m') 









Figure 5.2: The superblock consist of two enlarged block connected to each other. The 
two sites in the middle are the sites added last to the respective blocks. In the case of a 
Heisenberg chain the two enlarged blocks are connected only by the exchange of these two 
sites. 

The next step in the DMRG method is the formation of the superblock Hamiltonian (Fig. 
5.2). The superblock consists of two enlarged blocks connected to each other. In Fig. 5.2 
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open boundary conditions (OBC) are applied to the superblock. This BC are the most 
widely used in DMRG for it yields the best results for a given computational effort. We will 
discuss PEG latter. 

The DMRG method focus on a single eigenstate of the superblock Hamiltonian (usually 
the ground state), called the target state, which is used to construct the density matrix^. 
The ground state of the superblock is calculated (using Lanczos or any other method) . We 
then eliminate the states from the basis of the enlarged block that contribute the least to 
the ground state of the superblock. To calculate those, the density matrix is used. 

The concept of the density matrix was developed in statistical mechanics [26] by consid- 
ering the problem of a system in contact with a much larger bath. The ground state of the 
universe, i.e. system and bath, is known, and the question is which states of the system 
contribute the most to this ground state. This is what the density matrix can tell us. One 
can express the ground state of the universe (the superblock) in a basis that is the tensor 
product of the basis vectors of the system (one of the enlarged blocks) and the bath (the 
other enlarged block), 

mxD m'xD 

l*o)= E E a^M)®\b';)- (5.2) 

i=l j=l 

Hence many of the eigenstates of the system contribute to the one ground state of the 
universe. The density matrix of the system is given by 

m'xD 

Pw = E (5-3) 

We show an actual example of such a calculation below. The density matrix has the same 
dimension and block-diagonal structure of the Hamiltonian H^, for the enlarged block. If 
we denote by \ua) {a = 1, . . . ,m x D) the eigenstates of p and by Wa its eigenvalues then 
Y^aWa = 1 and Wa is the probability of the system being in the state \ua) given that the 
universe is in the state |^^'o)• 

This is the information we need to decide which states to keep in a renormalization group 
approach. In order to make an optimal decision of which states to discard and which to keep, 
it is a good criterion to consider the weight Wa of the states in the ground state of a larger 
system, which wc eventually want to describe. We must order the \ua) by their eigenvalues 
in a decreasing order and use the m of those states with largest eigenvalues to form a new 
basis for the enlarged block B^, which will then become B{1 + l,m). In symbols, 

HBii+i,m)^OH,0\ (5.4) 

where the rows of the m x {m x D) matrix O are formed by the \ua) previously selected. 
The change of basis in Eq. (5.4) renormalizes the Hilbert space, cutting its size back to m. 
Gonstructed in this way, the blocks are being prepared to be connected to another block in 
the next step, when a new superblock will be formed. By using the density matrix states 
we somehow look into the future and adapt the block for it. Besides Hb we also need to 

^It is possible to target several eigenstates simultaneously but, for a given computational effort, the 
accuracy decreases rapidly. 
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store other operators representing the sites at the border of the block. These operators 
are necessary to construct the interaction among the block and the site, when forming the 
enlarged block and also need to be transformed to according to Eq. (5.4). 

To illustrate the DMRG steps of enlargement and truncation, a full DMRG step for the 
antiferromagnetic spin 1/2 Heisenberg chain will be performed [2 7]. The starting point is a 
block 5(1, 2) of a single site. The possible states of the single site are 



IM = IT), l^'2) = |i). 



(5.5) 



For convenience the up/down basis is chosen. The basis itself is not stored. The only data 
that is stored are the operators needed to progress the algorithm namely, the operators 
needed to build the Hamiltonians for the enlarged block and the superblock. 

For one isolated site without external fields the Hamiltonian is zero. Since the up/down 
basis was chosen, the other operators are the spinmatrices given by 



(5.6) 



To build the enlarged system, another site is added. In this case the basis of the block is the 
same as the basis of the added site. 



|rfl) = |T), |^^2) = |i), 



(5.7) 



and the operators look the same as those of the block. Thus, the basis of the enlarged block 
is, by Eq. (5.1), 





= ITT) 


hi) 


= ITi) 


bt) 


= liT) 


K) 


= lU) 



(5.8) 



The Hamiltonian Hg for the enlarge block S(2,4) has non-zero elements, and describes the 

interactions of the sites in 5(2,4). He consists of the Hb, describing the interactions within 
the block, and the interactions between the rightmost spin of the block and the new site. In 
the above basis the Heisenberg Hamiltonian of the enlarged block is 



He = HB(Ejl 




(5.9) 
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and it looks as follows: 
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(5.10) 



In Eq. (5.9) the index h and d refers to the operators acting on the Hilbert space of the block 
and the site, respectively, and / is the unit matrix. In this first step we have m — D — 2 
but, as the block grows in size in the following steps we will have m > D. Note that only 
representations for the Hamiltonian of the block and for the operators S*"*", 5", and of 
the rightmost site of the block and the new site are needed to construct the enlarged block. 

The superblock is constructed by taking the enlarged block as the left block and connect- 
ing it to another enlarged block on the right (Fig. 5.2). In the so-called infinite size method, 
discussed in the next Section, the right block is the same as the left block, only spatially 
reflected so that the site last added to the left block is connected with the site added last to 
the right block. The rightmost site of the left block becomes the leftmost of the right block. 

In addition to the Hamiltonians of the enlarged blocks, one needs representations of the 
spin operators of the rightmost site of the enlarged block. In order to construct representa- 
tions for 5"+, S" , and in the basis of the enlarged block we have to calculate the tensor 
product of the unit matrix of the block Hilbert space and the operator in the representation 
of the basis of the rightmost site. For instance, the {S^)f, matrix, the ^'"'"-operator of the 
spin on the rightmost site in the basis of the enlarged block, is given by 




(5.11) 



Representations for {S^ )e and {S^)e are obtained in a similar way. The basis for the su- 
perblock is the tensor product of the bases from the two enlarged blocks being connected: 
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(5.12) 
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In general, and \h^) are distinct basis. Assuming that we want to calculate the ground 
state properties, it is possible to exploit 5"^ conservation and the fact that the ground state 
belongs to the subspace with 5"^ = 0. Therefore, we can concentrate only on states in this 
symmetry sector: 



(5.13) 
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The Hamiltonian of the superblock consists of three parts: the internal Hamiltonians of the 
two enlarged blocks and the exchange arising from the spin interacting at the connection 
between them: 



(s:), (s,-); + (5-). (s+);i + (53. » (53;, (s.m) 



where the prime refers to operators in the Hilbert space of the second enlarged block forming 
the superblock. We now build a representation for Hg in the basis l&f^^^): 
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(5.15) 



The ground state energy of is Eq — (1/4) (3 + 2-\/3) and the correspondent eigenvector 
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(5.16) 



In order to decide which states of the left enlarged block are the most important for the 
ground state of the superblock one uses the density matrix given by Eq. (5.3). Combining 
Eqs. (5.12) and (5.13) we obtain 
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(5.17) 
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This allows us to identify the coefficients Qij in Eq. (5.2) (and in Eq. (5.3)) and they are all 
zero except for: a 14, 022, 023, O32, 033, 041. For the density matrix we get 



12(2 + VS) 



/ 1 







ll + 6v^ 
-2(5 + 3^3) 
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(5.18) 



Note that p and (Eq. (5.10)) share the same block- diagonal structure. The eigenvalues of 
p are (1/12) (2 + ^3) ^ 0.02 for each of the triplet states and (21 + 12%/3) (12(2 + ^3)) ^ 0.94 
for the singlet state. The basis states are then ordered according to the size of the respective 
eigenvalues, with the singlet state (largest eigenvalue) coming first. The transformation 
matrix O in Eq. (5.4) is given by 



O 





(5.19) 



After determining the basis and the transformation, the representations of all operators used 
to describe the enlarged block are changed to the new basis. Applying the transformation 
to the He (Eq. (5.10)) leads to 



B{2,2) 



1-3 
4 I 



(5.20) 



In the present simple case and p have the same eigenvectors. That is the reason why 
the above transformation diagonalizes Hg. The same transformation is done with the other 
operators that will be need for future calculations. One example is the yS'"'"-operator, which 
has the following representation in the new basis 



(5.21) 



Note that, even though a site has been add to block 5(1,2) to form block B{2,2), the 
dimension of Hilbert space did not change due to the truncation performed. The states kept 
in the truncation are those with higher probability to be found in the ground states of the 
superblock system. 

In this example we have performed the truncation in the Hilbert space in order to illus- 
trate the procedure. In a practical calculation the system in our example would be too small 
to already start the truncation. Usually, we know from the beginning how many states will 
be kept. Thus, in the first steps the block is grown (sites are added) without truncation until 
the number of states needed to describe the block becomes larger than the number of states 
we want to keep^. If we had, for instance, decide to keep m = 20 states in the computation, 
the chain would be grown to a size of 5 sites with 2^ = 32 states, which is the first block 

''in this initial stops with no truncation the matrix O has dimension (m x D) x (m x D) and Eq. (5.4) 
becomes a simple change of basis. 
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size with Hilbert space dimension exceeding 20. Then, a truncation would create a block 
B{5, 20) from a enlarged block with 5 sites and in the following steps all blocks would have 
dimension 20, even though they represented a increasing number of sites. 

We often use the sum of the density matrix eigenvalues of the discarded states (1 — 
J2a=i ^ measure for the severity of the truncation. The goal is to keep this number as 

small as possible. In many cases it has been found that this number is roughly proportional 
to the error in the energy[28]. The proportionahty factor is of course model dependent. In 
doped fermionic models, we need to keep more states to achieve a good accuracy than in 
a spin model. Even for a given model the accuracy for a given truncation may depend on 
the parameters being used in the calculation {e.g., couplings and symmetry sector). For 
instance, close to a phase transition line or inside a critical (massless) phase the strong 
quantum fluctuactions tend to reduce accuracy. In the example above we discarded two of 
the triplet states leading to a truncation error of 0.04, which is unacceptably high. Truncation 
errors in actual calculations are usually kept smaller than 10^*^. 

In this Section the focus was on one DMRG step. Enlargement of the block by adding 
a site, the formation and diagonalization of the superblock, the calculation of the density 
matrix, and the truncation procedure were discussed. In the next Sections we describe how 
several DMRG steps are combined to calculate the properties of a given model. 

5.3 The Infinite System Algorithm 

The first implementation of the DMRG method was the infinite system algorithm[2, 3]. 
The goal was to use DMRG's advantage to decouple the system size and the size of the 
Hilbert space and calculate ground state energies of large systems, i.e., system size that are 
unreachable for exact diagonalizations, eventually converging to the thermodynamic limit. 

In the infinite system algorithm, the left enlarged block is connected by its own mirror 
image on the right side, so that the number of sites of the superblock is increased by two 
in each step. Growing the block and truncating the Hilbert space is done as explained in 
Section 5.2. Schematically the algorithm can be described as follows: 

1. Grow the chain to a size in which its Hilbert space dimension is just larger than m, 
the number of states to be kept. This is the first enlarged block. 

2. Form the superblock by adding an identical enlarged block on the right such that the 
sites, which were added last, are next to each other. 

3. Diagonalize the superblock, calculate the density matrix with respect to the left en- 
larged block. 

4. Change basis of the left enlarged block to the eigenbasis of the density matrix, keep only 
the m states with the largest density matrix eigenvalues. The transformed, truncated 
left enlarged block becomes the block for the next iteration. 

5. This block is enlarged, i.e., a site is added on the right. 

6. Continue with step 2 until convergence is reached. 
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Figure 5.3: Two DMRG steps of the infinite system algorithm, see text for discussion. 

In Fig. 5.3, two successive iterations of the infinite system algorithm are shown. The 
starting point is a block with / sites that is described by a basis with m states. Enlargement 
and construction of the superblock (step 2), leads to the situation portrayed in Fig. 5. 3. a. 

After diagonalizing the superblock, finding the transformation, and truncating of the 
enlarged block according to its density matrix (steps 3 and 4), one arrives at the situation 
depicted in Fig. 5.3.b. The new block describes a chain with I + 1 sites, but uses a basis 
with only m states. Enlarging the block (step 6) and building the superblock (step 2) leads 
to the situation in 5.3.C, and the procedure is repeated. 

From the computational point of view the most difficult part is the calculation and 
the subsequent diagonalization of the superblock Hamiltonian. The diagonalization can be 
done with the Lanczos method but any other method {e.g., the Davidson algorithm [2 9]) 
can be used. The computational effort depends on the size of the Hamiltonian matrix and 
accuracy needed for the ground state. Since the superblock Hamiltonian is block-diagonal, 
diagonalizing only the sector that has the proper quantum numbers reduces the matrix by 
a factor that depends on the of superblock, block, and model. Typically it is of order 10-20. 
Since the total Hilbert space for the superblock has dimension [DxmY, the most important 
determining factor for the size of the Hamiltonian is the number m of states kept in the block. 
In actual calculations m is typically a few hundred of states and is limited due to restrictions 
in computer memory and CPU-time. As we work with different models, the size D of the 
Hilbert space for a single site also affects the computational effort needed to reach a given 
accuracy in the results. 

To summarize: in the infinite system algorithm, the system size is increased in each step 
while the number of states kept to describe the blocks is constant. The goal is to grow the 
chain to a long-enough length, so that the energy and short range correlations around the 
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center have converged. The convergence is checked by keeping track of the difference AEq 
between the ground state energy of the super blocks in two sucessive steps^. 

5.4 The Finite System Algorithm 

In the finite system algorithm the goal is no longer to reach the thermodynamic limit, but 
rather to restrict ourselves to a finite system size L. In the beginning, until the superblock 
size reaches the system size, the algorithm is identical to the infinite size algorithm. When 
the system size reaches L, i.e. the enlarged block has L/2 sites, the left block is grown 
further but on the right a enlarged block with a smaller number of sites is used in order to 
keep constant the number of sites in the superblock. As soon as the decreasing size of the 
right block reaches a single site the procedure is stopped. One such iteration in which the 
left block has been calculated for all possible sizes, nearly up to L, is called a sweep over the 
system. After one sweep is done, we start all over with a small left block. However, from 
now on, the information about the best representation of the right block that complements 
the left block to the desired system size L is now present since it has been calculated in the 
previous sweep. When the optimal basis for a specific size of the left block is determined with 
DMRG, the result is stored and used in the next sweep as best possible guess for the optimal 
states describing the right block. The steps of the finite system algorithm are compiled in 
the following list: 

1. In the first sweep use the infinite size algorithm until the superblock size reaches the 
chain size L under investigation. After every truncation save all operators of the 
reduced block to disk. 

2. Enlarge the left block size I + 1. 

3. Read a block of size L — I — 2 from disk, this is the right block. 

4. Enlarge the right block to the size L — I — 1. 

5. Form the superblock from right and left enlarged blocks. 

6. Diagonalize the superblock, calculate the density matrix with respect to the left en- 
larged block. 

7. Change the basis of the left enlarged block to the eigenbasis of the density matrix, keep 
only the m states with the largest density matrix eigenvalues, save the block with the 
basis to disk. The transformed, truncated left enlarged block becomes the left block 
for the next step. 

8. Continue with step 2 until the right block becomes a single site 

9. If the right block is a single site begin a new sweep over the system i.e. construct a 
superblock with a left enlarged block containing two sites. Continue with step 3 until 
convergence is reached. 

'^{AEo)/2 converges to the ground states energy per site in the thermodynamic Hmit. 
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To illustrate the progress of the algorithm, two general steps are portrayed in Fig. 5.4. 
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Figure 5.4: Two DMRG steps of the finite system algorithm, see text for dicussion. 

At first glance Fig. 5.4 looks like Fig. 5.3. In fact the algorithms are very similar. The 
left bock is enlarged, the superblock is constructed and diagonalized (Fig. 5. 4. a). Then 
the density matrix for the left enlarged block is constructed and the block is reduced (Fig. 
5.4.b). This reduced block is then again enlarged, the superblock is built up and so on. 
The difference between the two algorithms is that the right block in Fig. 5.3 has the same 
number of sites I as the left side. In the finite system algorithm the right enlarged block 
always complements the left one to the target size L and thus becomes shorter, while the 
later grows. 

As an example, let us consider a calculation for the Hciscnbcrg chain with L = 16 sites 
and a truncation oi m = 24 states. The following superblocks have to be formed: 

• Initial sweep: 

[B(1,2),B(1,2)] [B(2,4),B(2,4)] [B(3,8),B(3,8)] [B(4,16),B(4,16)] 
[B(5,24),B(5,24)] [B(6,24),B(6,24)] [B(7,24),B(7,24)] [B(8,24),B(6,24)] 
[B(9,24),B(5,24)] [B(10,24),B(4,16)] [B(11,24),B(3,8)] [B(12,24),B(2,4)] 

• Following sweeps: 

[B(1,2),B(13,24)] [B(2,4),B(12,24)] [B(3,8),B(11,24)] 
[B(4,16),B(10,24)] [B(5,24),B(9,24)] [B(6,24),B(8,24)] 
[B(7,24),B(7,24)] [B(8,24),B(6,24)] [B(9,24),B(5,24)] 
[B(10,24),B(4,16)] [B(11,24),B(3,8)] [B(12,24),B(2,4)] 
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If the ground state of the superblock is even under reflection symmetry, left and right 
side can be interchanged. That means that the reduced left blocks stored during the first 
half of a sweep (when the size of the left block is smaller than L/2) can be already used 
as right blocks in the second half of the same sweep. However, if the ground state of the 
superblock does not show reflection symmetry, the reduced blocks from the left and right 
side have to be constructed and stored independently. Therefore, in the absence of reflection 
symmetry stored space and CPU-time is doubled. 

The finite system algorithm is terminated when convergence is reached, i.e. when the 
energy in succeeding sweeps does not improve (decrease) any more. It happens, however, that 
the energy stays on a certain level for two or three sweeps only to further decrease afterwards. 
Therefore, one cannot just compare the energies of the last two sweeps performed. The 
convergence behaviour of the model should be taken in account. This behaviour can be 
investigated by looking at small systems, where calculations are inexpensive. The number 
of sweeps necessary for convergence depends strongly on the system size L, the number of 
states kept m, and the model itself. Fermionic models need more sweeps than spin models, 
specially when they are doped. Typically, spin models converge in less then 10 sweeps, even 
for fairly large chains (L = 100 or larger), while calculations for the doped t — J model with 
the same m usually need roughly twice as many sweeps, even for smaller system sizes. 

The first sweeps do not yield very accurate results, their purpose is to generate a good set 
of blocks of different sizes. Therefore, the first sweeps are normally done with a small number 
of states. When convergence is approached the number of states kept in the truncation can 
be increased in order to improve accuracy. This helps saving CPU-time, specially when L is 
large. 

5.5 Measurement of Observables 

The ground state energy of the superblock is determined every time the superblock is di- 
agonalized. The value is used to determine whether convergence is reached. It turns out 
that in the finite system algorithm the energy is the lowest when the two blocks forming the 
superblock have the same size. Therefore, this symmetric configuration is used to measure 
also all other observables in which one is interested. 

Unfortunately the values for the other observables such as the value of a certain spin 
or the spin correlation between spins on different sites, are not as easily obtainable as the 
energy 

This is caused by the change of basis that is performed in every step. Even if we start out 
with a basis where the demanded properties of the basis states are known, or could easily be 
calculated, this knowledge fades fast with the repeated linear combination of basis states of 
new basis systems. Of course we could keep track, for instance, of {S^) for each site in each 
state but the computational effort would be enormous. 

The way that is chosen is to carry out the measurement by acturally evaluating the 
operator in the ground state. The expectation value of on site i is calculated as 

{St) = i^olSn^o). (5.22) 
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This is the textbook formula, but the difficulty of applying it here is not visible from this 
equation. The problem is the basis. The ground state is expanded in a basis that evolved 

in every step of the algorithm and could not have been anticipated at the beginning. 

When the site i is added to the block, the representation of the S'^-operator in the basis of 
the enlarged block is known (in Eq. (5.11) the calculation is done for S^). But in general this 
it too early to measure (Sf), because the symmetric configuration was not yet reached. In 
order to still have the right representation for in the symmetric configuration, the matrix 
has to be updated and stored every time the basis changes. Updating means that the basis 
change has to be performed on S'f . If {Sf)'j denotes the representation of the S'^-operator 
on site i in the Hilbert space of the enlarged block with j sites {j > i) and Oj is the matrix 
that transforms and cuts the basis before adding site j + 1, the representation of Sf after 
the truncation is 

{St), = 0,{StYp]. (5.23) 

When another site, the (j + l)th, is added, the representation of the operator also has to be 
adjusted according to 

{St)%, = {S!), ^ (5.24) 

and truncated with Oj+i. Following this procedure we always have representation of the 
operator in the current basis. 

When it is time to perform measurements, i.e. when both blocks forming the superblock 
have the same size, we only have to find the representation of the operator in the Helberb 
space of the superblock. If, e.g., the operator is acting on a site inside the left block this 
means tensorizing it with the unit element acting on all remaining spaces, namely, the two 
central sites and right block. 

If i is small, i.e., the site is close to the end of the block, a lot of basis changes have to 
be performed before the measurement is carried out. Due to the truncations that go with 
this procedure, the accuracy is decreased. We expect a greater accuracy from observables 
on sites close to the middle of the chain. 

The issue is somewhat more complicated for nonlocal operators, e.g. spin correlations 
like Cs{i,j) = {S^Sj). In general one could just take the representations of the involved S^- 
operators and multiply them, when the symmetric configuration is reached. However, there 
is a more accurate way to proceed. As an example, wc consider a spin-spin correlation where 
j = i + 1 and the symmetric configuration is reached at L/2 = i + 2. The representation of 
the two operators at the chain size i+2 are 

= {O,+im{h®S')Ol)®Ia)Ol,0h (5.25) 

Therefore, one gets for the spin correlation 

iStSt^^t+2 = {Oi+iimh ® S^)Ol) ® /,)Ol+i(0,+i(/, ® S^)Ol,) ® (5.26) 

Another way of calculating the matrix is to multiply the two operators as soon as possible. 
In this case the operation can be done when the enlarged block size is i + 1, which gives 

Cs{i, i + = Oi+,{{{Oi{h ® S')Ol) Id){h S'))Ol,. (5.27) 
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Then, from now on, {Cs{i,i + l))j+i is tranformed as a whole. Its representation for the 
enlarged block with i + 2 sites is 

(C,(z, I + = I + l))i+i ® Id. (5.28) 

Comparing Eq. (5.26) with Eq. (5.28) the difference is only a a Ol^^Oj+i-factor between the 
two S'^-operators. Without the truncation of states it would have been the product of two 
unitary matrices and thus be a unit matrix. The two ways of calculating {Cs{i, i+^))i^2 would 
be equivalent. With the truncation, however, this is no longer true, as one can immediately 
see calculating O^O with the projector in Eq. (5.19). The factor 0|_,_iOi+i leads to a loss of 
accuracy in the matrix multiplication since instead of multiplying the matrices, only their 
projections are multiplied. This error becomes worse the further the sites i and j in Cs{i,j) 
are apart, because every separating site add another pair 0^0. 

For DMRG procedure this means that we have to use the latter approach, multiplying 
the operators as soon as possible. Prior to the calculation a hst of the observables that we 
are interested on measure has to be made. When growing the chain, the on-site operators 
are stored as soon as they arc generated and updated every time the basis is changed. The 
products of two-site operators arc formed and stored as soon as on has a representation for 
both operators, then they are updated also updated. In the case of operators that involve 
more single site operators, e.g. pairing operators, we have to proceed in the same way. 

If we are measuring correlations between sites located on different blocks we can not 
multiply them before the symmetric configuration is reached. This means that measurements 
of correlations across the center will always have larger error than correlations where both 
sites are located in the same block. 

Prom these explanations it has hopefully become clear that measuring observables in- 
volves additional storage and calculations. In order to save computation time the measure- 
ment process is started as late as possible in the calculations, after convergence has been 
reached. In the infinite size algorithm we restrict ourselves to the sites close to the middle 
of the chain, which were the last to be added. In the finite size algorithm measurements are 
not performed in the initial sweeps when the energy has not converged yet. 

It is generally difficult to estimate the error on the values for the observables other than 
the energy for which it has been established that the error is proportional to the truncation 
error as discussed before. Unfortunately there is no known method to calculate the error 
of the observables from any other measured quantity. However, by checking how stable the 
results are as we change the number m of states kept in the truncation, we can have a 
qualitative control of the error. The error on the energy is generally smaller than that of the 
observables because it is a quantity that averages over the sum of many terms. 

5.6 General Remarks 

When working with fermionic models such as the Hubbard and t — J models we have to 
implement the anticommutation of the fermionic operators Eqs. (2.2) and (2.3) [30]. 

In order to implement PBC in DMRG a specific superblock configuration is required[3]. 
DMRG performs poorly under PBC. Typically, if a given accuracy is obtained under PBC 
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with m states kept in the truncation, then states will be needed to achieve the same 
accuracy under PBC[3]. 

An important improvement to DMRG involves keeping track of the wave function from 
step to step and perform a transformation into the basis corresponding to the current su- 
perblock. Since a good initial guess speeds up the Lanczos or Davidson convergence, this 
saves time in the diagonalization of the superblock[31]. 

When DMRG procedure converges to a fixed point the superblock ground state can be 
simply written as a matrix-product form and also be rederived through a simple variational 
ansatz making no reference to the DMRG construction. These very interesting analytical 
results are obtained in Ref. [32] and give some insight into the mechanisms working behind 
DMRG algorithms. 
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Chapter 6 

Conclusions and Acknowledgments 



The field of numerical simulations can stand by itself as a third way of doing science and its 
interaction with experiment and theory is very fruitful. We believe that scientific research 
in nearly every area in physics can benefit from it. As the power of nowadays computers 
increases rapidly the relatively new field of numerical simulations gains more and more 
prominence. To keep up with the advances in the hardware new methods and algorithms 
are being developed and traditional ones are being improved. The two techniques examined 
here are typical examples of such methods and algorithms. 

The author is thankful to those (too many to name here!) who carefully read this 
manuscript helping to improve it in many ways. The author also acknowledges hospitality 
at the Institute de Fisica de Sao Carlos - USP and the financial support from Fundacao 
de Amparo a Pesquisa do Estado de Sao Paulo - FAPESP - and Conselho Nacional de 
Desenvolvimento Cientifico e Tecnologico - CNPq - Brasil. 
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